Growth anisotropy of the extracellular matrix shapes a developing organ

Final organ size and shape result from volume expansion by growth and shape changes by contractility. Complex morphologies can also arise from differences in growth rate between tissues. We address here how differential growth guides the morphogenesis of the growing Drosophila wing imaginal disc. We report that 3D morphology results from elastic deformation due to differential growth anisotropy between the epithelial cell layer and its enveloping extracellular matrix (ECM). While the tissue layer grows in plane, growth of the bottom ECM occurs in 3D and is reduced in magnitude, thereby causing geometric frustration and tissue bending. The elasticity, growth anisotropy and morphogenesis of the organ are fully captured by a mechanical bilayer model. Moreover, differential expression of the Matrix metalloproteinase MMP2 controls growth anisotropy of the ECM envelope. This study shows that the ECM is a controllable mechanical constraint whose intrinsic growth anisotropy directs tissue morphogenesis in a developing organ.

Quantifications of pouch area in control discs (b) and in discs with reduced PPE growth (c) show that peripodial expression of PI3K DN does not affect DP growth. (e) Quantification of the pouch area (inner Wg ring) covered by mRFP positive peripodial cells. While in control discs the peripodial epithelium covers ~95% of the pouch, peripodial PI3K DP overexpression significantly reduces pouch coverage to ~66%. (f-g) Peripodial projection of control (e) and PI3K DN expressing discs (f) stained for E-cadherin (Ecad) to visualize cell outlines. Magnifications of the central PPE region marked by orange rectangles are shown top right (Z1 and Z2). (h) Quantification of peripodial apical cell area in control disc (black) and in disc where PPE growth was reduced by PI3K DN (red) as shown in (f-g). In control discs peripodial cell area is maximal in the centre of the disc and decreases towards the periphery. In PI3K DN discs the central PPE cell area is significantly reduced. (Error bands indicate standard deviation). Statistics: Statistical significance was assessed by a two-sided Student's t-test (unequal variance, *p≤0.05, **p≤0.005, ***p≤0.0005). In box plots the median is indicated by a central thick line, while the interquartile range (containing 50% of the data points) is outlined by a box. Whiskers indicate the minimum and maximum data range. Source data are provided as a Source Data file.
Y27632+Collagenase treated discs). n indicates number of discs. Statistics: Statistical significance was assessed by a two-sided Student's t-test (unequal variance, *p≤0.05, **p≤0.005, ***p≤0.0005). In box plots the median is indicated by a central thick line, while the interquartile range (containing 50% of the data points) is outlined by a box. Whiskers indicate the minimum and maximum data range. Source data are provided as a Source Data file.

Supplementary Figure 6 -Genetic degradation of the ECM results in epithelial thickness relaxation
(a) Two representative late 3 rd instar wing discs (disc 1 and disc 2) overexpressing Matrix-Metalloproteinase 2 (MMP2) in the posterior compartment (P-comp., magenta) for 12 hours before dissection in plane (left) and cross-section view (right). Temporal induction of MMP2 expression was controlled by a temperature shift to 29C thereby inhibiting a temperature sensitive version of the Gal80 (Gal80ts) protein. Discs were stained for Vkg (green) to visualize the ECM and F-Actin (Phalloidin) to mark epithelial outlines. While the posterior overexpression of MMP2 results in a loss of the ECM in the P-comp. and epithelial thickness relaxation, the anterior compartment (A-comp.) maintains its ECM and normal epithelial shape serving as internal control. (a') Magnification of the DP epithelium in the anterior versus the posterior compartment as indicated in (a, right: zoom 1 and zoom 2). (b) Quantification of epithelial thickness in the A-comp. (control) versus the P-comp. (no ECM). Genetic digestion of the ECM by MMP2 overexpression results in a reduction of epithelial thickness to ~20μm, a value very close to the reference thickness (~21.5μm) we observed upon acute digestion of the ECM using Collagenase (see Fig.3 and Supplementary Figure 5d). Therefore, the DP epithelium does not grow actively in z-direction but increase in tissue thickness are due to ECM mediated cell crowding and compression. Statistics: Statistical significance was assessed by a two-sided Student's t-test (unequal variance, *p≤0.05, **p≤0.005, ***p≤0.0005). In box plots the median is indicated by a central thick line, while the interquartile range (containing 50% of the data points) is outlined by a box. Whiskers indicate the minimum and maximum data range. Source data are provided as a Source Data file.  )). An increase in Vkg::GFP PPE peak density is observed from 72 to 118hAEL. Error band indicates the standard deviation. (b'') However, when ECM PPE thickness is quantified based on an intensity threshold (see dashed blue line in (b'), no change in thickness is observed. See methods for details on intensity profiles and thickness quantifications. (c) Representative section of the bottom ECM DP underlining the basal surface of the DP epithelium. A clear increase in thickness is visible from 72 to 118hAEL. This tendency to increase in thickness can be quantified in Vkg::GFP DP apicalbasal intensity profiles (c'). A significant increase in thickness is quantified in (c''). Error band indicates the standard deviation. (d) Plot of relative DP volume increase versus relative estimated ECM volume increase (obtained by multiplying the area of the inner Wg ring with the ECM thickness obtained in (c'')). We observe a linear correlation with a slope of ~0.8 (dashed line). This result suggests that the DP epithelium outgrows its EMC by ~20%. Error bars indicate standard deviation. Statistics: Statistical significance was assessed by a twosided Student's t-test (unequal variance, *p≤0.05, **p≤0.005, ***p≤0.0005). In box plots the median is indicated by a central thick line, while the interquartile range (containing 50% of the data points) is outlined by a box. Whiskers indicate the minimum and maximum data range. Relative increase in integrated Vkg::GFP PPE intensity versus volume increase of the peripodial epithelium (blue) and the DP epithelium (magenta, as quantified in Supplementary Fig.1k). Linear regressions are indicated by dashed lines, error crosses indicate standard deviation. Vkg::GFP PPE levels increase faster than the peripodial volume (slope of linear regression ~2.54) suggesting that the top ECM outgrows the peripodial cell layer. In contrast, the integrated Vkg::GFP PPE intensity and the DP volume increase to a similar extent (slope of linear regression ~0.95) from 72 to 118hAEL. (disc numbers PPE: n 72 =11, n 96 =10, n 118 =11, disc numbers DP: n 72 =9, n 96 =11, n 118 =12). (c) top: Maximum projections of the DP side of Vkg::GFP discs (green), stained for Wg/Ptc (magenta) at indicated stages. middle/bottom: Masked and extracted Vkg::GFP DP signal within the inner Wg ring underlining the DP. Projections in-plane (middle) and cross-section view (bottom). (c') Integrated fluorescent intensity of the Vkg::GFP DP underlining the DP (as shown in (c)). (d) Relative volume increase of the disc proper epithelium plotted versus the relative increase in integrated Vkg::GFP DP intensity (linear correlation, slope of ~0.85). This implies that the DP outgrows the ECM DP layer (also see Supplementary Fig.7d). (Linear regression is marked by dashed line, n 72 =9, n 96 =11, n 118 =12, error cross indicates standard deviation) (e) Assuming that the ratio of pouch volume increase divided by the ECM DP volume increase between 65h and 118h is 0.79, and that the growth of both layers is planar, we find that in our simulations this growth mismatch is not sufficient to describe the correct tissue geometry. Source data are provided as a Source Data file.

Supplementary Figure 9 -Changes in ECM morphology during larval growth (decellularization)
(a+b) ECM changes upon decellularization in fixed samples at different developmental stages in the PPE (a) and the DP (b). Fluorescent intensities and profiles are comparable between different timepoints and conditions. (a) top-middle: Cross-sections of the peripodial ECM labelled by Vkg::GFP in control (top) and decellularized (middle, Triton X-100 treated) wing discs at indicated age classes. Bottom: Profiles of peripodial Vkg::GFP levels (ECM PPE ) in control (black) and decellularized conditions (green) Error bands indicate standard deviation. (a') Quantification of peripodial ECM thickness at indicated time-points in control (black) and decellularized discs (green). During development we do not observe a significant change in ECM thickness, however, upon decellularization ECM thickness decreases. n indicates number of discs. (b) Same as in (a) but for the bottom ECM DP . (b') In contrast to the peripodial ECM, the thickness for the bottom ECM DP increases significantly during development. (c) Bottom ECM DP relative circular area plotted over time after the addition of Triton X-100 in ex vivo culture. The relative area decreases for 30min before reaching a plateau value at ~79% of the original area. (c') left: Cross-section of representative sections of the bottom ECM DP before and 60min after Triton X-100 addition. right: Quantification of Vkg::GFP DP intensity before and after decellularization. The profile of the relaxed bottom ECM DP shows increased peak intensity and increased width at the threshold value (dotted blue line, 10a.u.) chosen to quantify ECM thickness in Fig.4d''. (c'') Quantification of changes in estimated circular ECM volume (area A * thickness T) before (Vol 0 ) and 60min after decellularization (Vol 60 , n=26 circles / 9 discs) (d) Results for ex vivo decellularization in the ECM PPE layer. left: Three circular regions of interest (ROIs) were marked by photobleaching onto the ECM PPE of 118hAEL Vkg::GFP wing discs. right: Relative area changes after addition of Triton X-100. In contrast to the bottom ECM DP , the relative area in the top ECM PPE transiently increases before reaching a plateau at ~97% of original area after 60-70min. (d') ECM PPE thickness and Vkg::GFP density does not significantly change upon decellularization. (d'') Estimated ECM PPE volume marked by Vkg::GFP remains constant upon decellularization (n=20 circles / 7 discs). Statistics: Statistical significance was assessed by a two-sided Student's t-test (unequal variance, *p≤0.05, **p≤0.005, ***p≤0.0005). In box plots the median is indicated by a central thick line, while the interquartile range (containing 50% of the data points) is outlined by a box. . bottom: Average MMP2 profiles of control (black) and MMP2 KD discs (red). The domain between x=0 to x=-5 corresponds to the PPE (error bands indicate standard deviation). (d) Peripodial ECM posterior of the A/P boundary before (top) and after decellularization (bottom). Vkg::GFP PPE intensity profiles indicate ECM thickness (threshold value 50 a.u.) and peak intensity (maximum intensity of profile). (e) Relative change of ECM PPE thickness (Thickness 60min / Thickness 0min ) upon decellularization versus the distance between the peripodial and disc proper A/P boundary. While in control wing discs this distance is typically ~120µm at 118hAEL, in MMP2 KD this distance decreases with increasing knock-down efficiency. Consistently, the discs showing the strongest ECM thickness changes also show strong reduction in A/P compartment boundary distance. (f) Relative Vkg::GFP peak intensity changes (Intensity 60min / Intensity 0min ) upon decellularization plotted against A/P boundary distance. (g) Section view of control wing discs expressing RFP (green) in the posterior compartment. right, Magnification, peripodial height above the A/P DP is indicated by an orange arrow. (h) MMP2 KD wing disc expressing RFP and Mmp2 TRiP in the P-compartment. (i) Height of peripodial epithelium above the A/P DP . Statistical significance was assessed by a two-sided Student's t-test (unequal variance). In box plots the median is indicated by a central thick line, while the interquartile range (containing 50% of the data points) is outlined by a box. Whiskers indicate the minimum and maximum data range. (j-k) Late 3 rd instar wing disc expressing either CD8::RFP (magenta) alone (j, control) or together with a MMP2 TRiP line (k, MMP2 KD ) in the DP wing pouch (nub-Gal4), stained for Vkg (green). Central ECM DP is magnified to the right. (l) Quantification of basal disc outline. MMP2 KD in the DP results in slightly increased curvature. (m) DP knock-down of MMP2 leads to a significant increase in Vkg peak intensity and increased ECM thickness. These observations are consistent with the notion that MMP2 is also required in the DP tissue to fine-tune ECM growth anisotropy and disc shape.  M1 Description of the multi-layer wing disc model

M1.1 Theoretical framework
To gain insight into the shape of the growing wing disc and estimate the growth-induced stress, we modelled the Drosophila wing disc as a growing elastic tissue, implementing an existing theory of tissue growth that has been applied successfully to arteries and brain tissue [9,13,12]. The numerical implementation was solved in the finite element methods framework COMSOL Multiphysics. The tissue grows according to a specified growth deformation tensor, G. The total deformation gradient F consists of two contributions: The growth tensor G and the elastic deformation gradient A, i.e. F = AG. The multiplicative decomposition is explained in more detail in the caption of Fig. M2. We also made the assumption that the tissue is a hyperelastic nearly-incompressible neo-Hookean material [12,7,8] with the elastic energy W , given by [1,9,6] where µ and κ are the shear and bulk modulus of the material, respectively, and |A| is the determinant of the elastic deformation gradient A. Further,Ĩ 1 = I 1 |A| −2/3 , where I 1 is the first invariant of the right Cauchy-Green deformation tensor, I 1 = tr(A T A). Finally, the Cauchy stress tensor is given by To obtain the deformation of the body with a prescribed growth tensor G, we solve the balance of linear momentum div T = 0. The linear momentum balance is coupled with the morphoelastic decomposition F = AG, in conjunction with the constitutive law (M2). We assume that the external boundary of the wing disc is traction-free, Tn = 0, where n is the outward facing surface normal. The wing disc is assumed to be axisymmetric and we enforce a no-displacement boundary condition on one point along the symmetry axis. The full set of boundary conditions is shown in Fig. M1. In a cylindrical basis {E R , E θ , E Z }, we assume that the growth tensor takes the form where γ Z,i represents growth in the direction of the long cylinder axis and γ i growth in the direction perpendicular to it, i.e. in the circular plane which is the cross-section of the cylinder. If γ i = γ Z,i , * These authors contributed equally: Stefan Harmansa, Alexander Erlich † Correspondence: thomas.lecuit@univ-amu.fr no displacement stress-free axisymmetry   Figure M2: The wing disc is modelled as a multi-layer structure comprised of three layers or fewer. In the initial configuration B 0 , the tissue is ungrown and unstressed. The growth tensor G describes growth without stress, leading to an unstressed incompatible post-grown configuration B R , in which the individual layers have growth but do not fit into Euclidean space without breaking the connection between layers. Finally, the elastic deformation gradient A restores compatibility by introducing residual (internal) stress, bringing the body into the current (observed) configuration B t . The initial geometry in configuration B 0 are three glued-together discs of radius R 0 . The peripodial epithelium (PPE) is modelled as a disc of heigth H PPE , the disc proper (DP) as a disc of heigth H DP , and the extracellular matrix (ECM) as a disc with heigth H ECM . In the post-grown configuration B R , the three discs have radius γ Z, i R 0 and heigth γ Z, i H i with i ∈ {PPE, DP, ECM}, where the growth tensor in the cylindrical basis is G = diag (γ i , γ i , γ Z, i ).
growth is isotropic, and if γ i ̸ = γ Z, i , growth is anisotropic. Notice that we do not consider here the case of polar anisotropy, which means that we have implicitly assumed γ R,i = γ θ,i = γ i in (M3). We denote [G] the components of G in the cylindrical basis: The volume of the wing disc can be computed from the post-grown stress-free configuration B R , see Fig.  M2. This is done by adding the volumes of the three flat discs in B R : The volume of the deformed configuration B t will be the same in the limit κ → ∞, which corresponds to an incompressible material (|A| = 1).

M1.2.1 Scenario: Bilayer PPE-DP
In this scenario, we test the hypothesis that the wing disc consists of two layers, a fast growing peripodial epithelium and a slower growing disc proper epithelium. The ECM is not modelled in this scenario, see Table M1 column "PPE vs DP". For the sake of simplicity, we consider planar growth in both PPE and DP, that is The results of simulations with γ PPE = 4.3 and γ DP = 3.31, are shown in Fig. 3B. the relative volume increase of the PPE compared to the DP layer, normalized by initial volumes, is |G PPE |/|G DP | = γ 2 PPE /γ 2 DP = 1.69, in other words the PPE layer grew 70% more in volume than the DP layer. However, as discussed in the main text, experiments show that growth between PPE and DP is nearly compatible, demonstrate that epithelial doming and thickening are not due to a non-uniformity of growth within or between epithelial layers.

M1.2.2 Scenario: Bilayer DP-ECM, with differential growth anisotropy
In the main text, we describe how a two layer system made up of a DP layer and ECM layer capures the bending of the wing disc, as well as a number other of experimentally measured geometric quantities, see Fig. 5. In particular, the DP layer grows in plane whereas the ECM layer deviates from planar growth.
The growth tensors describing the two layers are The components of the growth tensor for the DP satisfy γ DP = 1 at t = 65h and γ DP = γ * DP at t = 118h and for the ECM they satisfy γ ECM = 1 at t = 65h and γ ECM = γ * ECM at t = 118h. For the disc proper, the parameter γ * DP can be determined by fitting the linear function which satisfies the above constraints. Since the volume over time in the DP, V DP (t), is available from experiments, we now show how γ DP (t) can be related to experimental volume measurements. The volume of the DP is given by V DP (t) = V 0,DP |G DP |, where |·| = det (·) denotes the determinant of a tensor. Here, |G DP | = γ 2 DP (t). The last two equations can be combined to The initial volume V 0,DP at t = 65h can be computed as the volume of a cylinder, V 0,DP = πR 2 0 H 0,DP where the values for the initial radius, R 0 , and initial heigth of the DP, H 0,DP are obtained from experimental measurements, as stated in Table M1 column "best fit DP vs ECM" (see the initial state, B 0 , in Fig. M2). With this information, we can solve (M9) for γ DP (t) = V DP (t) /V 0,DP . We make a least squares fit, using the experimentally measured values V DP (t) and V 0,DP and using the form of (M8) for γ DP (t). As a result of the least squares fit, we obtain γ * DP = 8.09.  Figure M3: Obtaining a relationship between γ * ECM and ρ, using volumetric data of the ECM. In the in-plane scenario ρ = 0, all of the volume increase in the ECM is assumed to be distributed in the plane, and the in-plane growth component γ * ECM is highest in this case. As ρ increases, an increased amount of the measured volume is assumed to be distributed in Z-direction, so that the in-plane growth component γ * ECM decreases. The special case ρ = 1 represents isotropic growth, in which case the ECM growth tensor takes the form G ECM = γ ECM 1, where 1 is the 3-dimensional identity.
To obtain γ * ECM , we proceed similarly. We use a linear functional form for γ ECM (t), which is identical to (M8) (with the subscript DP replaced by ECM). The ECM volume, on the other hand, depends on ρ, since |G ECM | = γ 2+ρ ECM (t). Thus, the ECM volume is given by where V 0,ECM = πR 2 0 H 0,ECM and values for R 0 , H 0,ECM are stated in Table M1 column "best fit DP vs ECM". With this information, we can solve (M10) for γ ECM (t) = (V ECM (t) /V 0,ECM ) 1 2+ρ . So for a given value of ρ, we make a least squares fit from experimental values V ECM (t) and V 0,ECM and using the linear form for γ ECM (t) given above. We repeat the least squares fitting for a series of values of ρ, obtaining data points (ρ, γ * ECM ) from the best fits, see Fig. M3. This data is well described by the smooth function γ * ECM (ρ) = 2.56 + 4.42e −1.33ρ .
With the parameters γ * DP , γ * ECM determined from volume data, the remaining parameters ρ and µ = µ DP /µ ECM are determined in Fig. 5B. There, we can see three distinct (partly overlapping) regions: In the blue region, labeled tol(H ECM ) ≤ 22.4%, the relative error between the mean experimentally measured value of ECM thickness and the simulation result for the corresponding pair of {ρ, µ ECM /µ DP } values does not exceed 22.4%. Similarly, the rose region labeled tol(H DP ) ≤ 22.4% compares in the same way experimental vs. simulated disc proper thickness. Finally, the orange region labeled tol(∆H/H ECM ) ≤ 3.53% considers the relative thickness increase after vs before decellularization, that is the ratio H after decell ECM /H before decell ECM . The same ratio can be relatively easily considered in the model, where ratio of reference ECM thickness to observed ECM thickness is γ ρ ECM H 0,ECM /H ECM (see Fig. 5a "reference state" and "observed state"). So in summary, in the orange region, the relative error between the mean experimentally measured value H after decell ECM /H before decell ECM and the simulation result γ ρ ECM H 0,ECM /H ECM for the corresponding pair of {ρ, µ ECM /µ DP } values does not exceed 3.53%.
In Fig. 5b, we can see that these three regions overlap in the dark shaded region, where the disc proper thickness, ECM thickness and ECM relative thickness increase upon decellularization are all simultaneously within the tolerance values when experimental and simulation results are compared. This region allows us to determine the best fit ρ = 0.45 and µ ECM /µ DP = 25.

M1.2.3 Scenario: Bilayer DP-ECM, both layers growing in-plane
This scenario serves as a demonstration that differential growth anisotropy is indeed essential to capture the wing disc morphology. In Supplementary Fig. 8e, using the parameter values given in Table M1 column "in-plane DP vs ECM", we explore the scenario from Section M1.2.2 but with a crucial difference: There is no differential growth anisotropy (ρ = 0), meaning that both the DP and ECM grow in-plane. The result is shown in In Supplementary Fig. 8e, demonstrating clearly that without growth anisotropy, the correct wing disc morphology can not be achieved, we end up with a structure that is far too flat.

M2 Numerical implementation
Numerical simulations are obtained using a finite element code that solves the equation of finite elasticity on a multi-layer structure of glued-together axisymmetric cylinders. The 2D axisymmetric numerical problem is solved in Comsol Multiphysics ® [3]. In Comsol, the morphoelatic decomposition F = AG cannot be directly entered into the software. We give here a succinct summary of the implementation procedure stated in works of Larry Taber [7,12], showing how growth problems can be studied in Comsol.

M2.1 General case
Comsol computes derivatives with respect to displacement gradients to obtain a second Piola-Kirchhoff stress tensor where we denote C F = F T F the right Cauchy-Green strain tensor of the total deformation gradient F. In particular, the finite element formulation requires the stress S per unit initial area (B 0 ). The second Piola-Kirchhoff stress in terms of Cauchy stress is given by For an incompressible or nearly incompressible material, the Cauchy stress is given by (M2), which can be rewritten as where we denote C A = A T A the right Cauchy-Green strain tensor of the elastic deformation gradient A. Inserting (M13) into (M12), we get As shown in [7], the gradient of W satisfies the following transformation relationship: Inserting this into (M14), we find So the appropriate expression of the second Piola-Kirchhoff stress can be obtained by multiplying the equation for S, which Comsol uses by default, by |F| / |A|. In addition, W must be defined in Comsol in terms of the components of C A , that is .

M2.2 Specific case
Note that Comsol will solve for C F . Therefore, it is necessary to provide Comsol with the components of C A in terms of the components of C F , in order for a successful implementation of the modified strain energy W and the modification (M16). In our problem, we assume that the growth tensor is written in a cylindrical basis {E R , E θ , E Z } where it has a diagonal form PPE vs DP (Fig. 2) best fit DP vs ECM (Fig. 5 These can be inserted into the re-defined strain energy density for Comsol, (M17), and used to compute the determinant |A| for the re-defined second Piola Kirchhoff stress (M16).

M3 Validation of the numerical implementation against exact solution
In this section, we test the computational framework presented in Section M2 by comparing it with an analytically known solution for an incompressible neo-Hookean material that grows. The setup we are considering is shown in Fig. M4A. We consider a cylinder with an initial radius R 0 , with no displacement allowed in Z-direction, that is [F] ZZ = 1. In the analytical axisymmstric calculation, we impose [F] ZZ = 1 everywhere in the bulk, whereas in the numerical implementation, as shown in Fig. M4A, we impose this constraint on the top and bottom highlighted planes. Further, we impose a polar growth anisotropy. Our goal is to compare the numerical and analytical expressions for the Cauchy stress tensor and to determine the relative error between them, so that we establish a benchmark for what kind of error to expect in the simulations presented in Section M1.2 where no analytical solution exists.

M3.1 Derivation of exact solution
We consider the case of a single incompressible growing neo-Hookean disc. We assume that there is no deformation at the point of symmetry, and that there are no external forces, so that any deformation is caused purely by growth and the elastic response. This calculation follows a similar path to [5,2,6]. For an incompressible disc in the cylindrical basis {E R , E θ , E Z }, the morphoelastic decomposition F = AG (M23) Eliminating α, we obtain the kinematic relationship r (R) r ′ (R) = γ R γ θ R, r (0) = 0 .
Let W (α R , α θ ) be the strain-energy density, which relates to the Cauchy stress tensor by where p is the Lagrange multiplier enforcing incompressibility. In components this reads where we used the shorthand notation [T] RR = T R , [T] θθ = T θ . With no external loads, mechanical equilibrium requires div T = 0, which takes the form ∂T R /∂r = (T θ − T R ) /r. Expressing with respect to the reference radius R, this becomes Defining W (α) := W α −1 , α , we have Now taking into account (M28), for the radial stress we must solve For a neo-Hookean strain-energy density and taking into account the relationships from the morphoelastic decomposition (M2) and (M24), we can express (M29) as Once the radial stress is known, the circumferential stress can be obtained from (M28) To obtain the full solution, one must provide γ R (R) and γ θ (R) and can then solve the two non-linear coupled ODEs (M28) and (M31). For γ R , γ θ spatially constant, we obtain the exact analytical solution Figure M4: Comparison between analytical solution and numerical implementation in Comsol.

M3.2 Comparison of exact and analytical solutions
In Fig. M4, we compare the exact analytical solution (M34), (M35) with the numerical implementation in Comsol Multiphysics ® as described in Section M2. We obtain excellent agreement, with the largest relative error between numerics and analytics at the disc center being 3.36% and the larges error at the disc edge being 1.38%. Given that at the disc center, the stress tensor has a singularity (T ∼ log R) as see in (M34) and (M35), we expect the lower error at the edge to be more representative. The Finite Element mesh in this case was composed of 26006 triangular elements (all simulations presented in Section M1.2 had around 25k triangular elements, always using the Comsol predefined setting "Extremely fine" for the mesh, without any manually defined mesh refinement). Overall, the comparison between exact and analytical solutions shows that the numerical implementation is reliable, and justifies its use for more complex geometries like the sandwich structures presented in the paper.